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Abstract 

We have derived exact Langevin equations for a model of quasispecies dynamics. The 
inherent multiplicative reaction noise is complex and its statistical properties are 
specified completely. The numerical simulation of the complex Langevin equations is 
carried out using the Cholesky decomposition for the noise covariance matrix. This 
internal noise, which is due to diffusion-limited reactions, produces unavoidable 
spatio-temporal density fluctuations about the mean field value. In two dimensions, 
this noise strictly vanishes only in the perfectly mixed limit, a situation difficult to 
attain in practice. 



1 Introduction 



The study of replicator models of prebiotic evolution has received considerable 
theoretical attention during the past three decades [1] . There is a renewed and 
active interest in this subject owing to the crucial fact that the dynamics of 
real viral populations is known to be described by quasispecies [2,3,4]. To date, 
the bulk of the theoretical work devoted to the analysis of quasispecies dy- 
namics has tacitly assumed spatially homogeneous conditions [1,2,3,4,5]. The 
importance of diffusive forces has begun to be recognized and taken into ac- 
count [6,7], but rather less attention, if any, has been payed to the presence of 
the unavoidable internal density fluctuations [8] that are necessarily present in 
all realistic incompletely mixed and diffusing systems of reacting agents [9] . In 
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the deterministic approach to chemical and molecular species that diffuse and 
react, the fluctuations are simply ignored. Nevertheless, it is well known that 
if the spatial dimensionality d of the system is smaller than a certain upper 
critical dimension d c , the intrinsic fluctuations do play a decisive role in the 
asymptotic late time behavior of decay rates, and the results obtained from 
the mean field equations are not correct [10]. Even far from asymptotia, these 
fluctuations can strongly affect the dynamics on local spatial and temporal 
scales [11]. The mean field limit is strictly valid only in the infinite diffusion 
limit, because the reactions themselves induce local microscopic density fluc- 
tuations that can be amplified by the underlying nonlinear dynamics. This 
internal noise is multiplicative, and scales as the square root of the product 
of concentrations, with the consequence that the noise does not necessarily 
vanish in the large particle- number limit. 

Given the relevance of molecular replicator models such as the quasispecies for 
viral dynamics and the key role played by fluctuations in nonlinear systems, 
this Letter has a twofold objective. First, we outline a rigorous derivation of the 
exact stochastic partial differential equations (SPDE) that govern the dynam- 
ics of a simple quasispecies model (denotes in this work a system comprised of 
a master plus mutant molecular species in competition for limited resources). 
Second, we study the nature of the internal fluctuations on the replicator evo- 
lution. But in order to achieve this latter goal, we must also develop and apply 
a method suitable for simulating multi-component SPDE's with complex noise 
[12]. The validity of this approach is confirmed by numerical simulation. Nei- 
ther of these objectives has been achieved before. We illustrate this method 
with a simple model that is analytically tractable, but the technique given here 
can be applied to more realistic chemical and molecular reaction processes. 



2 Model 



We consider a molecular quasispecies model with error introduced via the 
faulty self-replication of the master copy into a mutant species. The mutant 
species or, error-tail, undergoes non-catalyzed self-reproduction, but has no 
effect on the master species. The system is closed, only energy can be ex- 
changed with the surroundings, where activated monomers react to build up 
self-replicative units. These energy rich monomers are regenerated from the 
by-product of the reactions by means of a recycle mechanism (driven by an 
external source of energy) maintaining the system out of equilibrium. The 
closure of the system directly imposes a selection pressure on the population. 
In what follows, M*,I,I e denote the concentrations of the activated energy 
rich monomers, the master and the mutant copies, respectively. The kinetic 
constants are introduced in the following reaction scheme. 
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Accurate replication with rate AQ: 

M * + jAQ 2L ^ 



Error replication with rate A(l — Q): 

M* + I A( ^3 ] I + I e . (2) 



Mutant-species replication with rate A e : 

M * + I e ^ 2I e . (3) 



Degradation of the master and mutant copies into monomers with rates r,r e , 
respectively: 

I M, I e M. (4) 



Re-activation of energy-depleted monomers (induced by a generic energy- 
driven reaction): 

M e ^M*. (5) 



The quality factor Q e [0, 1]. In order to keep the following development 
mathematically manageable, we will assume that the monomer reactivation 
step Eq.(5) proceeds sufficiently rapidly, so that we can in effect, regard the 
decay of / and J e , Eq.(4), plus the subsequent reactivation M — > M* as 
occurring in one single step, i.e., 

IAM^M^I^M* (6) 
le J^ M er^y M *^i e A+ M *. (7) 

Here, r', r' e are correspondingly, the combined effective decay plus regeneration 
rates of I,I e to the reactivated monomers M*, although in what follows, we 
will simply write r and r e . This means that our continuum field theory, will 
depend on three (M*, J, I e ), instead of four (M, M*, I, I e ) concentration field 
variables. If we suppose the system is being bathed continuously by an external 
energy source, the monomer reactivation Eq. (5) is occurring continuously, and 
this should be a reasonable approximation. To complete the specification of 
the model, we allow for spatial diffusion. This is incorporated into the master 
equation associated with the above reaction scheme. The diffusion constants 
are denoted by -D S ,D/ and D e for M, I and I e , respectively. The constraint 
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of constant total particle number is automatically satisfied by the continuous 
chemical concentration fields in the mean-field limit. Most importantly, this 
provides a selection pressure on the quasispecies. 



3 Complex Langevin equations 

It is straightforward to write down the continuous time master equation corre- 
sponding to the above reaction scheme Eqs. (1,2,3,6,7). This master equation is 
then mapped to a second-quantized description following Doi [13]. From Doi's 
operator language, we then pass to a path integral representation in terms of 
continuous stochastic fields [14] to obtain the following action S valid for any 
space dimension d (D. Hochberg, M.-P. Zorzano and F. Moran, unpublished 
notes): 



In the mean field limit the continuous field variables a,b,c in Eq.(8) cor- 
respond to the physical densities of the molecules M, J, I e , respectively. The 
fields a, b, c, are conjugate to the internal fluctuations. Only if this noise is real, 
do a, b, c continue to represent the physical densities. Otherwise, for imaginary 
or complex noise, these fields will likewise be imaginary or complex, but their 
stochastic averages (a), (b) , and (c) do however correspond to the real physical 
densities [15]. 

Since the a, b, c fields appear quadratically in S, we can make use of the 
Hubbard-Stratanovich transformation to integrate over them exactly in the 
path integral / VaVaVbVbVcVce^ s[a ~ a ^ c ' t] . This final step yields a prod- 
uct of delta-functional constraints which directly imply a set of exact coupled 
stochastic partial differential equations with specific noise properties (see be- 
low). 

In particular, in a two-dimensional space and making the reasonable assump- 
tion that both master and mutant species diffuse equally D Ie — Dj — D 
and have equal effective degradation plus reactivation rates, i.e. r e /r = 1, the 
reaction-diffusion system in dimensionless variables 1 is given as follows: 



Define the dimensionless fields: a = —a, b = —b, and c = — c; the dimensionless 




D s V 2 a + Aab + A e ac — rb — r e c) 
AQab + rb) 

A(l — Q)ab — A e ac + r e c) 

A(l — Q)ab(bc — ab) — A e ac(c 2 — ac 



(8) 
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= (D s /D)V 2 a -ab-ac + b + -c + rj a 
Or o 

— = V 2 b + Qab-b + fj b 
or 



dc 

where V 2 



V 2 c + 5(1 — Q)ab + Sac — c + r) ci 



+ and rj = (r] b , f) c , r) a ) is the noise with (17) 



and 



correlation matrix B = (ffff ): 



B(x,y; t) 



eQab ^e<5( 1 — Q)ab — \eab 



\ 



\e5(l - Q)ab e5 2 ac 



V 



— \eab 



— \eb~ac 



— \e5ac 
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The initial condition a , b , c and the ratio of replication rates 5 = A e j A < 1 
defines the scaled number of particles in the closed system N — J J dxdy(a + 
fro + tF)- In the two-dimensional case the total number of particles is given by 
the ratio N = N/e where e = AjDj is the ratio of the reaction to the diffusion 
processes (in any dimension d, we have e = (r/Di) d / 2 A/r). In the absence of 
noise ff this react ion- diffusion system has a set of homogeneous solutions: (i) 
b = c = 0, a — N, (the trivial solution); (ii) 6 = 0, c/S — N — k, a — |, if 

5 > 1/N and 6 < Q; and (iii) c/5 = ( %^ Q \ b = ^^g"^ , a = 1/Q 
if Q > 5 and Q > 1/N. We want to understand how this deterministic, 
mean-field homogeneous solution, is modified by the internal multiplicative 
noise term with defined spatio-temporally varying covariance (the fields vary 
in both space and time). From inspection of the covariance matrix B we see 
that this noise is proportional to e, the ratio between the reaction and diffusion 
processes. In the limit e — > (infinite diffusion, i.e., a perfectly homogeneous 
system) there are no fluctuations, and we recover the mean field result. The 
noise term vanishes when b(x, t) = c(x, t) = (the trivial solution is an inactive 
state) and/or when a(x,t) = 0. If the system stays close to the mean field 
result, then a — 1/Q and b and c scale as (NQ — 1). Therefore the noise 
covariance is expected to scale with eN. Our main interest here is the limit 
when the stochastic term makes a significant contribution eN = e 2 N > 1. In 
this regime the problem can not be analyzed by perturbation theory and thus 
must be treated numerically 2 . 



( \ 1/2 

time r = rt and spatial coordinates £j = ( -/jj 1 Xj. Their corresponding derivative 
operators are ^ = ^t|, and V 2 = (^r)'V 2 . The dimensionless noises are defined by 



_ _ ±4, and V 2 , 

or r at ' v r 

Va = ^Va,i)b = ^Vb, and f) c = jfrj c . 
2 After the space and time discretization of the stochastic partial differential equa- 
tions, the numerical integration of the finite-difference equations has been performed 
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Note that B is a symmetric matrix with det£> = — ja 3 bce 3 5 2 (b + c), i.e. it is 
negative definite. Thus it either has one or three negative eigenvalues, suggest- 
ing that at least one noise component will have negative auto-correlation, and 
will therefore be imaginary. We expect this situation to be relatively common 
in reaction-diffusion systems with multiple species. 



4 Results 



We next apply the Cholesky decomposition, which is used when the symmetric 
matrix is positive definite, to relate this correlated noise to an uncorrelated 
one. We apply this algorithm for the first time here to with negative 

definite correlation matrix to obtain the "square root" L, where some of the 
terms are manifestly imaginary. This is very useful for relating the noise rf to 
a new real white Gaussian noise £ with = 1, such that rf = L£ 3 : 



- ( 5 (1- Q)y/b 6y/4Qc-(l-Q)*b 

+ 1 (l-Q)b-2Qc 
2 VQ 1 2 VQy/4Qc-{l-Q) 2 b 2 

+ Vea _ £ 3 . (9) 

UQc-{l-Q)% 




This transformation allows us to separate the real and imaginary parts of the 
noise, a extremely useful feature to have for setting up and running a stable 
numerical simulation. 

The noise rf a for the nutrient field a, Eq.(9), always has an imaginary com- 
ponent and since all the equations are coupled to the nutrient a, even if the 
initial condition for a, b, c is real, the fields will, in principle, have both real 
and imaginary parts, and thus we have to solve a system of 6 partial differen- 
tial equations, one for each field (JRa, £ra, £sb, 3?c, Sc) with two-dimensional 
diffusion and noise. We expect the imaginary fields to be zero in the average, 



using forward Euler with time step Ar = 5 x 10" 4 and a spatial mesh of 154 x 154 
cells with cell size Ax = Ay = 0.35 and periodic boundary conditions. 
3 For an M x M symmetric matrix B, we can apply the Cholesky Decomposition 
B = LL T to extract the square root of the matrix in the form of a lower triangular 



matrix L with L u = B u - Y?k=i L fk and L j* = T~i ( B ji ~ Sfc=i L ik L jk) for j = i + 
1,...,M. Note that (ffff T ) = (L$ T L T ) = L(££ fT )L T = LL T = B is automatically 
satisfied. 
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since the stochastic averages (a), (b) , and (c) correspond to the physical den- 
sities [15]. 

Integrating the system numerically we confirmed that the scaled number of 
particles was conserved and remained real KiV = N, and SsN = (within 
computational errors) 4 . We found that, independently of the initial condi- 
tion and noise intensity, it always converged to the same state with spatial 
fluctuations around this value: (^) x ,y — ^ V< ^~(i-I)~ < ^ , (b) x ,y — ^qTI^T 1 ^ ■> 
{a)x,y — 1-1 Q, which is the mean-field solution 5 An example of these spatial 
fluctuations is shown in Fig. (1) at r = 25 for e = 1 and N = N = 10 3 , 
with D s /D = 10, Q = 0.92, and 5 = 0.8. In the case of b and c, the densities 
show spatial fluctuations of up to +1% (in white) and —1% (in black) with 
respect to the mean field value. Inspection of the figure clearly shows that 
the master copy and mutant species are spatially anti-correlated. In the case 
of a, the fluctuations range in the order of ±5% about the mean field value 
a — l/Q. Next we explore the dependence of the spatio-temporal field fluctu- 



a(x,y,T — 25) b(x, y,r = 25) c(x,y,r = 25) 




Fig. 1. Long-time spatial distribution of the fields when 
Q = 0.92,5 = 0.8, D s /D = 10 and e = 1, N = N = 10 3 . The three fields 
show fluctuations greater than 1% about the mean-field value. 

ations cr%,, = ((&(t) — (&(T))x,j/) 2 )x,# etc., in e and N. In Fig. 2-left we show 
the deviation scaled by y/Ne for iV = ^ = 10 3 and two different cases of 
e. As expected, ai scales with %/iVe, in this 0.04 x y/~Ne. Next we show 

the same for N = 10 4 , the mean value of erg is the same but the deviations 
with respect to it (which is related to the kurtosis of the distribution) increase 
with N . Notice that the relative weight of these fluctuations with respect to 

4 In fact, for computational purposes, QiV was set as small as possible. We con- 
firmed that both KiV = / / dxdy($ta + m + ^a) = N and = f f dxdy(%a + 
96 + ^) = Nx 10~ 7 are conserved during the integrating time and that the results 
were independent of 9?iV. 

5 The imaginary part converged to {Qc/S) X; y = ^pzr~^ ; {^b)x,y = "^ffg ^ ; 
(3>a)x,y = 0. The spatial average is denoted by (•)£,§, which by the ergodic hy- 
pothesis, is equivalent to averaging over the noise. 
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the mean-field value decreases with the number of particles since b oc N and 
therefore a~ b /b oc 0My/Ne/(Ne) oc O.OA/VN. 




Fig. 2. Time evolution of the fluctuations u^r) = V ((^( r ) ~~ (K r ))) 2 ) scaled by 
y/Ne = y/We for N = 10 3 (left) and N = 10 4 (right) with e between 1CT 3 and 1. 
The fluctuations around the mean- field value increase with N and e. 



5 Conclusions 



Deriving the Langevin equation description from the Master equation provides 
a quick analytic survey of the global system dynamics (stationary states, etc), 
as well as the strength and characteristics of the internal noise and its de- 
pendence on the spatial dimensionality. For any dimension, the system can 
be solved numerically with high accuracy. This is in contrast to the direct 
Monte Carlo simulation of the system, such as in [16], where no preliminary 
analytic information can be given on the effects of reaction noise or dimen- 
sionality dependence. We have shown that the local deviations with respect 
to the mean-field solution scale with the number of particles and the ratio be- 
tween the reaction and diffusion processes. The fluctuations vanish only in the 
perfectly mixed limit, when the diffusion processes are infinitely fast (e = 0). 
Thus, even in the case of very high diffusion rates, an autocatalytic system can 
show important deviations with respect to the perfect mixing limit, provided 
that the number of particles in the reactor system is sufficiently high. For sys- 
tems with a higher degree of non-linearity (i.e., third-order), these deviations 
may eventually lead the system to new asymptotic states and also induce the 
formation of spatial patterns [8]. Regarding the error catastrophe, this concept 
differs from virology to molecular replicator dynamics. In the former, the error 
catastrophe implies the proliferation of lethal mutations with the subsequent 
loss of viral infectivity, and consequently, the extinction of the viral population. 
This is reflected in our model as the trivial solution (i). On the other hand, 
solution (ii) corresponds to the error catastrophe state predicted by molecular 
replicator dynamics where there is no wild type population, while there is an 
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error tail. We have found that this solution is an absorbing barrier. Finally, 
spatial diffusion and internal noise must in fact play an important role in the 
correct description of viral infection dynamics, where typically, bottlenecks of 
different intensities may lead to infection of host cells by limited numbers of 
viral particles [17]. 

We thank Carlos Escudero for many useful discussions and for working through 
some preliminary analytic calculations and Esteban Domingo for reading the 
manuscript and suggesting additional references. M.-P.Z. is supported by an 
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